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Abstract 

A  3D  simulation  tool  for  modeling  solid  oxide  fuel  cells  is  described.  The  tool  combines  the  versatility  and  efficiency  of  a  commercial  finite 
element  analysis  code,  MARC®,  with  an  in-house  developed  robust  and  flexible  electrochemical  (EC)  module.  Based  upon  characteristic 
parameters  obtained  experimentally  and  assigned  by  the  user,  the  EC  module  calculates  the  current  density  distribution,  heat  generation, 
and  fuel  and  oxidant  species  concentration,  taking  the  temperature  profile  provided  by  MARC®  and  operating  conditions  such  as  the  fuel 
and  oxidant  flow  rate  and  the  total  stack  output  voltage  or  current  as  the  input.  MARC®  performs  flow  and  thermal  analyses  based  on 
the  initial  and  boundary  thermal  and  flow  conditions  and  the  heat  generation  calculated  by  the  EC  module.  The  main  coupling  between 
MARC®  and  EC  is  for  MARC®  to  supply  the  temperature  field  to  EC  and  for  EC  to  give  the  heat  generation  profile  to  MARC®.  The  loosely 
coupled,  iterative  scheme  is  advantageous  in  terms  of  memory  requirement,  numerical  stability  and  computational  efficiency.  The  coupling 
is  iterated  to  self-consistency  for  a  steady-state  solution.  Sample  results  for  steady  states  as  well  as  the  startup  process  for  stacks  with 
different  flow  designs  are  presented  to  illustrate  the  modeling  capability  and  numerical  performance  characteristic  of  the  simulation  tool. 
©  2004  Elsevier  B.V.  All  rights  reserved. 
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1.  Introduction 

A  solid  oxide  fuel  cell  (SOFC)  is  a  device  that  converts 
chemical  energy  into  electrical  energy  continuously  by 
means  of  an  electrochemical  (EC)  process.  The  high  oper¬ 
ating  temperature  of  SOFCs  has  numerous  advantages.  For 
example,  kinetic  activity  is  high  and  nickel  is  a  sufficiently 
good  catalyst,  eliminating  the  need  to  use  costly  precious 
metals.  CO  can  be  used  as  fuel  instead  of  as  a  poison,  as 
in  other  low  temperature  fuel  cells,  offering  potential  for 
fuel  flexibility  and  high  efficiency.  Natural  gas  fuel  can  be 
reformed  within  the  cell  stack,  eliminating  the  need  for  an 
expensive  external  reformer  system.  SOFCs  also  offer  the 
possibility  of  heat  co-generation  for  even  higher  efficiency. 
However,  the  high  operating  temperature  can  lead  to  com¬ 
plex  materials  problems,  including  mechanical  stress  due  to 
the  different  thermal  expansion  coefficients  of  the  cell  com¬ 
ponents;  electrode  sintering,  interfacial  diffusion  between 
electrolyte,  electrode  and  interconnect  materials.  All  these 
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factors  decrease  the  efficiency  and  the  stability  of  the  cell. 
The  choice  of  component  materials  is  limited,  and  fabrica¬ 
tion  problems  further  increase  the  cost  of  manufacturing. 
SOFCs  also  require  more  time  and  energy  for  startup.  Such 
problems  have  limited  the  commercial  development  of 
SOFCs. 

Planar-type  designs  have  received  much  attention  recently 
because  they  potentially  offer  higher  power  density  than  the 
tubular- type  SOFC.  This  is  attributed  to  the  low  electrical  re¬ 
sistance  due  to  shorter  current  paths.  In  particular,  a  very  thin 
electrolyte  film  can  be  used  in  an  anode- supported  SOFC, 
drastically  reducing  the  Ohmic  resistance  and  enabling  op¬ 
eration  at  intermediate  temperatures,  maintaining  this  bene¬ 
fit  while  alleviating  the  associated  material  and  mechanical 
difficulties. 

Despite  wide  interest  and  intensive  research  as  well  as 
significant  progress,  current  SOFCs  represent  an  immature 
technology.  To  be  commercially  feasible,  major  develop¬ 
ment  is  required  on:  (1)  lowering  the  cost  of  manufacture; 
(2)  improving  durability  to  match  that  of  other  electrical 
generation  techniques  in  use;  and  (3)  increasing  efficiency 
to  approach  the  theoretical  limit,  thus  broadening  the  range 
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of  potential  applications.  This  development  requires  exam¬ 
ining  a  huge  parameter  space  of  potential  operating  temper¬ 
atures,  materials,  plate  and  component  designs,  geometries, 
dimensions,  optimization  of  stack  configurations,  fuels,  re¬ 
formers,  heat  exchangers,  and  others. 

Although  the  design  and  operation  of  an  SOFC  appear 
deceptively  simple,  many  of  the  phenomena  dominating  the 
performance  of  the  SOFC  are  complex,  competing,  and  not 
well  understood.  Experimental  advances  are  limited  because 
physical  prototyping  is  expensive,  time  consuming,  and  can¬ 
not  explore  a  range  of  designs  and  operating  conditions.  To 
speed  up  the  development  and  understanding  of  SOFCs,  it 
is  highly  desirable  to  establish  mathematical  models  that  in¬ 
corporate  the  known  physics  and  behavior  of  SOFC  materi¬ 
als  to  predict  and  improve  performance. 

A  mathematical  simulation  of  an  SOFC  is  helpful  in 
examining  various  operating  parameters  such  as  tempera¬ 
tures,  materials,  geometries,  dimensions,  fuels,  reformers, 
etc.  and  determining  their  associated  performance  charac¬ 
teristics.  Such  a  simulation  can  be  used  to  determine  the  ef¬ 
fect  of  varying  the  design  and  operating  parameters  on  the 
power  generated,  fuel  conversion  efficiency,  maximum  cell 
temperature  reached,  stresses  caused  by  temperature  gradi¬ 
ents,  and  the  effects  of  different  coefficients  of  thermal  ex¬ 
pansion  for  the  electrolyte,  electrodes  and  interconnect.  In 
addition,  such  simulations  can  provide  answers  to  questions 
such  as  how  much  the  electrical  properties  of  the  cell  ma¬ 
terials  need  to  be  improved  or  what  the  air  and  fuel  flow 
rates  must  be  to  avoid  excessive  temperature  and/or  pressure 
drop.  Thus,  mathematical  simulation  offers  the  potential  to 
direct  the  technology  development,  test  the  significance  of 
various  design  features  and  effectiveness  of  developments 
in  materials  or  fabrication  procedures,  and  select  optimum 
operating  conditions  from  the  feasible  set  of  process  param¬ 
eters.  Modeling  of  SOFCs  is  a  critical  aspect  of  the  SOFC 
technology  development  process. 

A  key  feature  of  the  Core  Technology  program  spon¬ 
sored  by  the  US  Department  of  Energy  is  the  development 
of  modeling  and  simulation  tools  that  are  applicable  to  a 
wide  range  of  SOFC  designs  and  easily  accessible  to  the 
SOFC  research  and  development  (R&D)  community.  As 
the  constitutive  thermal,  chemical,  EC,  and  transport  pro¬ 
cesses  are  strongly  coupled,  the  SOFC  operations  involve 
multi-physics  processes  and  require  versatile  multi-physics 
tools  for  realistic  descriptions.  However,  no  available  com¬ 
mercial  code  has  the  capability  to  treat  the  fully  coupled 
physical  processes  involved  in  SOFCs.  Though  there  are 
many  model  development  efforts,  most  resulting  proprietary 
software  has  significant  limitations  and/or  not  widely  avail¬ 
able.  One  exception  is  the  software  developed  by  Recknagle 
et  al.  [1]  who  created  a  SOFC  modeling  tool  based  on 
a  commercial  computational  fluid  dynamics  (CFD)  pack¬ 
age,  STAR-CD.  However,  the  CFD  algorithm  of  STAR-CD 
makes  the  numerical  stability  and  computational  efficiency 
less  favorable.  To  assist  the  technology  development  that 
involves  design  optimization  of  various  geometric,  material 


and  operational  parameters,  parametric  studies  are  often 
required.  As  the  cost  of  such  parametric  studies  increases 
exponentially  with  the  number  of  working  parameters,  and 
many  such  parameters  are  involved  in  SOFC  designs,  the 
computational  efficiency  of  the  modeling  tool  is  of  critical 
importance  in  many  regards  and  is  highly  desirable.  To  this 
end,  a  numerically  stable  and  efficient  algorithm  is  preferred. 

As  commercially  available  finite  element  (FE)  codes  are 
well  developed  and  have  a  long  history  of  successful  appli¬ 
cation  in  structural  design,  commercial  FE  codes  appear  to 
be  natural  choices  for  the  platform  of  the  SOFC  modeling 
tool.  In  particular,  the  MARC®  FE  code  from  MSC  Soft¬ 
ware  is  a  multi-physics,  robust  structural  (mechanical,  ther¬ 
mal,  etc.)  analysis  tool.  MARC®  has  a  full  library  of  user 
subroutines  that  allow  for  user  interaction  with  the  code. 
Its  user-defined  functions  allow  algorithms  describing  the 
chemical  and  EC  processes  and  heat  generation  to  be  in¬ 
corporated  efficiently.  Consequently,  MARC®  was  selected 
as  the  underlying  platform  for  integration  with  the  in-house 
developed  EC  module,  which  models  the  EC,  chemical  re¬ 
action,  and  heat  generation. 

The  goal  of  this  paper  is  to  describe  the  capability  of 
our  EC  module  and  the  coupling  of  our  EC  model  with  the 
MARC®  FE  analysis  (FEA)  code.  The  EC  model  determines 
the  current  distribution,  gas  species  concentration,  and  heat 
generation,  and  the  FE  code  solves  the  heat  transfer  prob¬ 
lem.  This  tool  will  allow  for  a  variety  of  designs,  materi¬ 
als,  geometries,  and  flow  conditions  to  be  evaluated  at  both 
transient  and  steady-state  operating  conditions. 

Section  2  of  this  paper  describes  the  theoretical  back¬ 
ground  and  associated  technical  details  of  the  EC  software 
module.  The  EC  theory  contains  a  detailed  account  of  the 
chosen  current-voltage  (/-V)  relations  and  associated  heat 
production.  Determination  of  the  equilibrium  contents  of 
multi-component  fuels  is  also  described  to  treat  fuels  with 
CO  and  CH4  and  possibly  other  hydrocarbon  fuels.  The 
calculation  of  the  heat  production  due  to  the  chemical 
reactions  is  also  described.  Section  3  describes  the  over¬ 
all  algorithm  and  functionality  of  the  EC  module  and  the 
linkage  between  the  EC  module  and  the  MARC®  FEA 
code.  Section  4  shows  a  prototype  SOFC  cell  and  stack 
model  and  the  corresponding  sample  results.  The  numerical 
performance  of  the  coupled  MARC®-EC  package  is  also 
described.  A  brief  summary  is  given  in  Section  5. 

2.  Theoretical  background  and  technical  details 

In  a  fuel  cell  operation,  the  flow,  thermal,  chemical,  and 
EC  systems  are  intrinsically  coupled.  The  temperature  and 
flow  field  determine  the  rate  of  chemical  and  EC  reaction. 
Heat  generation  and  absorption  by  the  chemical  and  EC  re¬ 
actions  in  turn  affect  the  temperature  distribution  and  gas 
flow  composition.  The  coupling  requires  an  iterative  proce¬ 
dure  to  achieve  a  self-consistent  solution.  MARC®  performs 
the  thermal  analysis  if  all  heat  flux  is  properly  accounted 
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for.  The  key  linkage  between  the  EC  module  and  MARC® 
is  for  the  EC  module  to  feed  MARC®  with  the  heat  flux  re¬ 
sulting  from  the  chemical  and  EC  reactions,  which  are  based 
on  the  temperature  field  MARC®  provides  the  EC  module. 
The  description  for  the  heat  transfer  solution  is  referred  to 
in  the  MARC®  documentation,  and  only  the  basic  theory 
and  implementation  details  relevant  to  the  EC  module  are 
described  here. 

The  EC  process  generates  the  electrical  power  needed  for 
the  SOFC  applications.  SOFC  efficiency  depends  on  opti¬ 
mization  of  plate  design,  stack  configurations,  and  operating 
conditions  as  well  as  understanding  the  complex  interactions 
between  various  EC  phenomena.  The  EC  and  chemical  pro¬ 
cesses  generate  heat,  which  results  in  thermal  stresses  and 
affects  structural  reliability.  Consequently,  the  EC  process 
is  the  heart  of  the  fuel  cell.  One  of  the  most  important  as¬ 
pects  of  SOFC  design  is  the  voltage  and  current  distribution 
in  the  cell  plate.  This  couples  with  the  temperature  distribu¬ 
tion  from  the  thermal  and  flow  model  and  also  with  the  EC 
reactions  at  the  electrodes. 

Overall,  the  EC  module  gives  the  solutions  to  the  current 
density  distribution,  species  concentration,  chemical  reac¬ 
tion,  heat  generation,  and  flow  field.  The  EC  is  treated  at 
the  effective  property  or  continuum  level,  i.e.,  the  simula¬ 
tion  of  electrode  and  cell  performance  involves  some  pa¬ 
rameterized  EC  models.  Such  an  EC  model  is  described  as 
a  current-voltage  relation,  or  I-V  curve,  for  a  single  cell,  in 
terms  of  parameters  that  are  effective  cell  properties  and  op¬ 
erational  parameters.  The  solutions  to  the  chemical  reactions 
(CO-water  shift  reaction,  CH4  internal  reforming,  etc)  are 
based  on  the  equilibrium  theory.  The  flow  solution  is  based 
on  the  approximation  of  a  quasi-two-dimensional  laminar 
flow  (flow  in  each  gas  channel  is  treated  as  one-dimensional). 
Mass  balances  are  used  with  the  flow  pattern  to  establish 
species  concentrations  and  fluxes  at  any  point  in  the  fuel 
cell.  Each  aspect  of  the  EC  module  is  described  below. 

2.7.  I-V  model  for  H2  fuel 

In  the  EC  module  of  the  MARC®-EC  package,  the  EC 
is  based  on  continuum-level  I-V  relations.  The  I-V  rela¬ 
tion  describes  the  voltage  (potential)  loss  at  a  specified  cur¬ 
rent  with  respect  to  the  ideal  thermodynamic  performance, 
which  is  called  overpotential  or  polarization.  That  is,  the 
I-V  relation  is  a  model  for  calculating  current  density,  cell 
voltage,  and  heat  production  in  an  SOFC  cell  with  H2  or 
other  fuels,  taking  as  inputs  local  values  of  the  gas  partial 
pressures  and  temperatures.  This  cell  I-V  curve  is  specific 
for  the  materials,  structural  characteristics,  and  operational 
parameters  (gas  compositions,  pressure,  temperature)  of  a 
given  PEN  (cathode-electrolyte-anode)  element.  The  user 
can  change  the  relevant  parameters  to  reflect  their  cell’s  EC 
performance. 

Two  I-V  formulas  and  the  associated  parameters  are  im¬ 
plemented  in  the  EC  module.  In  this  section,  we  describe  in 
detail  the  I-V  relation  that  fits  best  with  H2  fuel,  though,  as 


discussed  later,  it  can  also  be  used  for  composite  fuel  without 
modification  in  most  cases.  Implementation  details  are  also 
discussed.  Generally,  the  I-V  relation  can  be  written  as  [2]: 


Here,  C  =  RT/4F ,  R  is  the  gas  constant,  T  the  tempera¬ 
ture,  and  F  (96487  Coulomb/mol)  the  Faraday  constant.  In 
Eq.  (la),  Eopen  is  the  equilibrium  (open  circuit)  voltage,  or 
emf  (electromotive  force),  of  the  cell,  i  the  current  density, 
iRi  the  Ohmic  potential  drop,  and  io  the  exchange  current 
density;  b  sinh_1(//2/o)  represents  the  electrode  activation 
polarization,  io2  the  O2  transport  limiting  current,  z'h2  the 
H2  transport  limiting  current,  p the  H2  partial  pressure, 

and  the  H2O  partial  pressure  at  the  fuel  channel.  The 
last  three  terms  represent,  respectively,  the  O2,  H2  and  H2O 
concentration  polarization.  Notice  that  b  sinh_1(//2/o)  re¬ 
duces  to  the  familiar  Tafel  form  for  i  z'o,  7>ln(///o),  while 
eliminating  its  logarithmic  divergence  at  i  =  0.  All  the 
parameters  in  Eq.  (la)  have  clear  physical  meanings  and 
their  dependencies  on  the  temperature,  cell  geometry,  and 
gas  partial  pressures  are  known  by  experimental  data  and 
transport  theory.  In  the  following,  the  theory  and  the  eval¬ 
uation  detail  of  each  quantity  that  is  necessary  for  software 
implementation  are  described. 

The  thermodynamic  cell  potential,  or  the  Nernst  potential, 
£open,  depends  on  reactant  and  product  partial  pressures  as 
well  as  temperature.  For  example,  for  fuel  with  hydrogen 
only, 


E 


open  — 


-AG 

IF 


AG  is  the  free  energy  change  of  the  reaction  H2  +  1  /2O2  — > 
H20: 

AG  =  G3  -  Gi  -  ±G2  (3) 


where  G  is  the  free  energy  per  mole  molecules.  Subscript 
1-3  refer  to  H2,  O2,  and  H2O,  respectively. 

The  pressure  and  temperature  dependence  of  the  Nernst 
potential  can  be  determined  according  to  the  thermodynam¬ 
ics.  As  ideal  gas  is  a  good  approximation  for  the  species  in¬ 
volved  in  the  operating  condition,  the  pressure  dependence 
of  the  Nernst  potential  can  be  written  as: 


A  Gp0  is  the  free  energy  change  when  all  of  the  species 
are  at  Po  =  1  atm  and  is,  therefore,  called  the  standard  cell 
potential  or  standard  emf.  It  depends  only  on  temperature. 
T\  and  T2  are  the  fuel  and  oxidant  temperature,  respectively. 
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Using  the  approximation  Cp  =  a  +  bT  -\-  cT 2,  which  is 
found  to  fit  the  experimental  data  of  many  gases  to  within 
0.5%  over  T  from  0  to  1250  °C  [3],  the  temperature  depen¬ 
dence  of  the  standard  free  energy  can  be  calculated  as: 


G(T,  Po)  =  G(T0,  Po)  +  a{T  -  T0)  +  ^(T2  -  T2) 

+  fT3-  T3)  -  [TS{T,  Po)  ~  T0S(T0,  P0)l 

(5) 


where 


S(T,Po)  =  S(T0,  P0)+alnf 

^0 


+  b(T 


To)  +  C-{T2 


and  a ,  b ,  c,  S(Tq,  Pq),  and  AG(To,Po)  are  experimental  data 

[3]. 

The  electrical  resistance,  R(,  of  the  PEN  structure  is  cal¬ 
culated  simply  as  the  sum  of  the  cathode,  Rp,  electrolyte, 
Re ,  anode,  Rn ,  and  bipolar,  R\>,  electrical  resistances: 

Ri  =  Rp  Re  Rn  +  Rb  (2) 


The  ionic  resistance  of  electrolyte  is  calculated  as  [4] : 


Re  =  l  qCqT  exp 


where  /e  is  the  thickness  of  the  electrolyte  and  Qq  the  ef¬ 
fective  activation  energy.  Rp,  Rn,  and  R\>  are  assumed  to  be 
Ohmic  and  have  linear  dependence  on  the  temperature.  The 
temperature  dependence  of  the  electrical  resistance  can  also 
be  defined  by  the  user  if  desired. 

For  the  thin  cell  structure  we  are  interested  in,  Rp+Rn+Rb 
is  much  smaller  than  RQ  and  may  be  neglected  in  prac¬ 
tice.  However,  these  terms  are  explicitly  calculated  in  the 
software  so  the  tool  can  describe  cells  with  thicker  elec¬ 
trodes/interconnects  or  other  materials  with  some  desirable 
properties  but  lower  conductivity. 

Notice  that  there  are  many  other  factors  of  practical  im¬ 
portance  not  considered  above.  These  include:  (a)  contact 
resistance:  this  depends  on  manufacturing  process  and  can 
very  well  be  more  than  the  above  resistance;  (b)  rib  dimen¬ 
sion:  as  channels  increase  the  resistance  of  interconnect;  (c) 
formation  of  highly  resistive  interfacial  products  such  as 
La2Zr207  and  SrZr03  during  cell  operation  and  the  elec¬ 
trode  adhesion  period  [5].  One  can  add  such  data,  if  desired, 
in  computing  parameters  such  as  resistance,  current  distri¬ 
bution,  and  heat  generation. 

Generally,  the  electrode  activation  polarization,  Z?sinh-1 
(//2/o),  is  strongly  temperature  dependent.  The  electrode  ac¬ 
tivation  polarizations  are  controlled  by  the  reaction  kinetics 
of  the  respective  electrodes.  They  represent  the  voltage  loss 
incurred  due  to  the  activation  necessary  for  charge  trans¬ 
fer  and  are  usually  dominated  by  the  cathode  process.  One 
implementation  in  the  package  is  to  describe  the  temper¬ 
ature  dependence  by  interpolation  and  extrapolation  based 
on  the  available  experimental  data  at  several  temperatures 


provided  in  [2].  Notice  that  these  data  are  cell  dependent, 
and  corresponding  experimental  data  should  be  used  in  real¬ 
istic  simulations.  Another  implementation  is  to  express  the 
activation  polarization  in  a  Butler- Volmer  form,  i.e.,  writ¬ 
ing  b  =  RT/aF ,  where  a  is  the  Butler- Volmer  transfer  co¬ 
efficient.  The  exchange  current  density  corresponds  to  the 
dynamic  electron  transfer  rate  at  equilibrium,  which  is  ther¬ 
mally  activated.  Therefore,  the  exchange  current  density  can 
be  expressed  as  [6]: 

/(,  =  Px  exp  (9) 


where  the  prefactor,  Px,  and  the  activation  energy,  Eact,  are 
properties  specific  for  the  electrode-electrolyte  interface  in 
question. 

Regarding  the  concentration  polarization,  the  H2  limiting 
current  density  is  calculated  as: 


2FDeff(T)p°H2 
(RTl  a) 


(10) 


where  Deff(7)  is  the  effective  H2  diffusion  coefficient,  ^ 
the  H2  partial  pressure,  and  la  the  anode  thickness.  In  one 
option,  the  temperature-dependent  portion  of  the  expression 
is  found  by  interpolation  and  extrapolation  of  the  experi¬ 
mental  data.  Another  option  is  to  determine  the  temperature 
dependence  according  to  the  transport  theory  (described  in 
more  detail  in  Section  2.2).  A  similar  method  is  used  to  find 
the  O2  limiting  current  density,  though  it  is  usually  incon¬ 
sequential  for  an  anode- supported  cell. 


2.2.  I-V  model  for  fuel  with  CO 


The  anode  concentration  polarization  terms  of  Eq.  (la) 
imply  that  all  current  is  produced  by  H2  oxidation.  For 
fuel  with  CO,  the  measured  current  is  the  sum  of  the  cur¬ 
rents  derived  by  H2  (z’i)  and  CO  (if)  oxidation.  Accordingly, 
the  current  for  the  anode  concentration  polarization  terms 
in  Eq.  (la)  should  be  replaced  with  i\ .  Similar  adjustment 
also  applies  to  the  current  generated  by  CO  oxidation.  More 
specifically,  Eq.  (la)  should  be  changed  to: 


V(i)  =  Eopen  -  iRi  -  b  sinh  1  l  )  +  c  In  /  1 


2/0 


10: 


,  1  -  i\  \  f  1  +  i\ 

+  2C  In  I  — — -  I  —  2C  In  7 


lU: 


E0pen  IRi  b  sinh 


*h2o 


-1 


2/0 


+  Cln 


1  -  i 

b2 


,  1  -  h  \  (\  +h 

+  2C  In  I  - ^|-2Cln/ 


ICO 


1  co2 


(lb) 


where  /h2o  =  Pu2o^  / Rib  l^e  Ibnftmg  current  for  H2O 

transport,  ico  and  ico2  =  Pcofco/Pco  are>  respectively, 
the  limiting  current  for  CO  and  CO2  transport;  b!  is  used  to 
indicate  the  possible  difference  on  the  effective  activation 
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polarization  due  to  the  CO  current  path.  Generally  speaking, 
CO  is  used  in  two  ways,  one  by  oxidation  and  produces 
current,  the  other  by  shift  reaction  and  produces  H2.  If  the 
current  is  produced  largely  by  CO  oxidation,  a  reverse  shift 
reaction  could  occur  and  H2  (and  CO2)  could  be  converted  to 
CO  (and  H2O).  On  the  other  hand,  if  the  current  is  produced 
largely  by  H2  oxidation,  shift  reaction  would  compensate 
for  the  depleted  H2.  Though  the  same  amount  of  current 
would  generate  a  similar  amount  of  heat  and  resulting  similar 
gas  composition  by  the  two  processes,  the  latter  scenario 
is  more  likely  as  shift  reaction  is  energetically  favorable 
(9.84kcal/mol).  In  fact,  it  is  estimated  that  about  98%  of 
current  is  produced  by  H2  oxidation  in  common  situations 
[7].  This  implies  that  Eq.  (la)  can  also  be  used  for  CO  fuel 
with  minimal  error.  That  is,  we  can  often  assume  current  is 
produced  only  by  H2  oxidation. 

CO  is  shifted  by  the  equilibrium  condition.  Parameters 
used  for  H2  fuel  can  also  be  used  here  without  modifica¬ 
tion  (though  one  should  bear  in  mind  that  it  would  be  very 
different  if  CO  is  responsible  for  most  of  the  current).  This 
simplification  is  adopted  in  the  software  as  one  option.  An¬ 
other  option  is  to  explicitly  consider  the  current  produced 
by  CO  oxidation.  The  total  current  by  the  simultaneous  H2 
and  CO  oxidation  is  treated  by  a  parallel  electrical  circuit 
analogy  [6].  Assuming  the  anode  activation  overpotential  is 
negligible  compared  to  the  cathode  activation  polarization, 
or  if  the  activation  overpotential  for  H2  and  CO  oxidation  is 
the  same,  and  the  shift  reaction  does  not  proceed  before  the 
oxidation  reaction,  the  diffusive  polarization  should  be  the 
same  for  H2-H2O  and  CO-CO2  transport: 


-77ac  =  2Cln 


=  2Cln 


—  2Cln 

-  2Cln 


(11) 


With  some  algebra,  the  solution  to  Eq.  (11)  with  a  given 
total  current  yields: 


t  =  exp 


-B+  (B2-4AC)1/2 
2 A 


(12) 


where 

i  1  1 

A  —  7 - 7 - +  7 - +  7 -  (13a) 

iu2oico2  *h2o  i  co2 


i  i  11 

B  —  7 - 7 —  +  7 — 7 - 7 —  +  7 — 

*h2o*co  iu2ico2  iu2  i  co 

i  1  1 

C  =  7 — : - t - — 

in2ico  in2  ico 

The  current  components  are  given  by: 
1  -t 

(1/ fit2)  +  it/  /h2o) 


1  1 

z’h2o  ico2 

(13b) 

(13c) 


(14a) 


1  -  t 

(1/  A:o)  +  it/  ico2 ) 


(14b) 


These  equations  give  the  I-V  relation  when  CO  is  present 
in  the  fuel  stream.  To  complete  the  description  for  the  limit¬ 
ing  current  as  indicated  by  Eq.  (10),  one  needs  to  determine 
the  effective  diffusion  coefficients  entering  into  the  limiting 
currents  for  each  gas  species  in  the  anode. 

Diffusion  coefficients  for  binary  gas  mixtures  are  fairly 
well  known  and  scale  fairly  well  with  the  result  of  elemen¬ 
tary  kinetic  theory: 

D  oc  p/2{M\~x  +  M2-1)1/2 

-  (15a) 

pa 


where  T1^2  is  the  reduced  mass  reflect  a  rms  thermal  speed, 
a  factor  Tip  comes  from  the  density,  and  a  a  scattering 
cross-section.  We  use  a  similar  formula  due  to  Fuller  et  al.  [8] 
that  incorporates  a  weak  T-dependence  of  the  cross-section: 


10— 7T1'75(jMi_  1  +M2_1)1/2 

pin  +  0')2 


(15b) 


where  D  is  in  m2/s,  T  in  Kelvin,  the  molecular  weights  are  in 
amu,  the  pressure  is  in  atmospheres,  and  the  r’s  are  tabulated 
phenomenological  molecular  radii  in  arbitrary  units.  Data 
from  [8]  are  used. 

The  full  theory  for  a  multinary  gas  mixture  involves  a  dif¬ 
fusion  matrix  with  off-diagonal  elements  in  which  the  gradi¬ 
ent  of  one  component  drives  part  of  the  current  of  the  others. 
In  addition  to  the  mathematical  complexity,  the  off-diagonal 
matrix  elements  are  not  well  known.  As  an  alternative,  we 
compute  a  “unitary”  diffusion  coefficient  for  each  species 
based  on  its  binary  coefficient  with  each  of  the  others.  As 
it  is  scattering  rates  or  inverse  conductivities  that  tend  to  be 
additive  in  such  situations,  a  reciprocal  average  weighted 
by  the  fractional  concentration  y*  of  each  gas  constituent  is 
used: 

(A-r1  =  y>,(A,r1  (i6) 

j 


The  binary  coefficients  in  Eq.  (15)  should  be  replaced 
with  the  unitary  coefficients  in  Eq.  (16)  for  the  limiting  cur¬ 
rent  calculations.  In  the  limit  of  an  arbitrary  binary  mixture, 
Eq.  (16)  gives  the  exact  result.  It  also  approaches  a  ternary 
mixture  correctly  when  one  constituent  is  very  dilute.  With 
these  formulae,  the  anode  concentration  polarization  can  be 
calculated  for  any  given  current  and  the  I-V  relation  is  de¬ 
termined. 


2.3.  Model  for  composite  fuel:  chemical  equilibrium 

Though  most  tests  done  for  SOFC  are  done  with  H2  fuel, 
fuel  with  CO,  CH4,  or  other  hydrocarbon  may  be  desirable 
in  practice.  The  ability  to  treat  fuels  other  than  H2  is  required 
for  modeling  flexibility.  The  effects  of  CO  fuel  on  the  I-V 
relation  and  the  possible  treatments  have  been  described 
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above  in  detail.  The  explicit  determination  of  the  composi¬ 
tion  change  of  the  fuel,  however,  is  still  required.  Moreover, 
more  discussions  are  needed  for  hydrocarbon  fuel.  The 
general  model  we  used  in  the  software  is  that  the  current  is 
either  purely  produced  by  H2  oxidation  (see  Section  2.1)  or 
by  both  H2  and  CO  oxidation  (see  Section  2.2).  The  hydro¬ 
carbon  is  used  not  by  oxidation  but  by  steam  reforming.  The 
assumption  is  valid  in  most  cases,  though  there  are  some 
exceptions  [9,10].  The  fuel  composition  change  is  deter¬ 
mined  by  assuming  chemical  equilibrium.  The  equilibrium 
theory  does  not  address  the  dynamic  issue,  i.e.,  how  fast  a 
system  achieves  equilibrium,  which  is  strongly  material  and 
process  dependent.  While  CO  shift  reaction  is  a  fast  process 
and  can  be  assumed  always  to  be  at  equilibrium,  the  steam 
reforming  reaction  may  be  a  slow  process  and  may  not  have 
reached  equilibrium  in  the  fuel  cell  operation  [11].  The  de¬ 
gree  of  deviation  from  the  equilibrium  and  the  effect  on  the 
modeling  results  depends  on  the  actual  fuel  cell  operating 
characteristics.  More  general  treatment  will  be  implemented 
later.  The  description  below  discusses  the  details  of  deter¬ 
mining  the  fuel  equilibrium  compositions  in  the  EC  module. 

2.3.1.  Fuel  with  CO 

For  a  given  initial  partial  pressure  of  the  CO  fuel,  Ph 2, 
Ph2o,  Pco,  and  Pco2,  the  determination  of  the  resulting 
content  requires  four  conditions.  Since  no  nuclear  reaction 
is  involved,  each  of  the  total  contents  of  H,  C  and  O  should 
remain  constant.  The  initial  condition  provides  three  con¬ 
straints  on  the  resulting  content: 

Ph  =  2(Pr2  +  Ph2o)  =  constant  (17a) 

Pc  =  Pco  +  Pco2  =  constant  (17b) 


As  indicated  by  the  simple  analytical  expression,  the  equi¬ 
librium  content  determination  requires  negligible  extra  com¬ 
putation  time. 

2.3.2.  Equilibrium  composition  for  fuel  with  CH4 

For  fuel  with  CH4,  there  are  five  undetermined  composi¬ 
tions,  H2,  CO,  H2O,  CO2,  and  CH4.  It  requires  one  more 
constraint  than  that  of  CO  fuel.  This  additional  constraint 
is  provided  by  the  equilibrium  constant  of  the  CH4  steam 
reforming  reaction: 

CH4  +  H20  ->  CO  +  3H2  (22) 

Unlike  the  case  with  CO  fuel,  it  is  difficult  to  derive  simple 
analytical  expressions  for  the  equilibrium  compositions  in 
this  case.  However,  it  is  easy  to  devise  a  numerical  solution. 
The  algorithm  for  determining  the  equilibrium  compositions 
involves  the  following  steps: 

(1)  Initial  estimate  CH4  content. 

(2)  Providing  a  sensible  H2,  CO,  H2O,  CO2  composition, 
obeying  the  total  content  of  H,  C  and  O  as  the  input 
and  get  the  resulting  shift  equilibrium  compositions,  as 
described  in  the  last  subsection. 

(3)  Check  to  see  whether  the  steam  reforming  equilibrium 
condition  is  satisfied  within  the  desired  accuracy.  If  not, 
go  back  to  Step  (1)  until  the  equilibrium  content  is 
found. 

The  software  implementation  for  these  steps  is  quite  sim¬ 
ple,  and  the  overall  equilibrium  content  can  be  obtained 
easily. 

2.4.  Flow  model 


Po  =  Pco  +  2Pco2  +  Pr2o  =  constant  (17c) 

Thus,  the  remaining  three  compositions  can  be  easily 
found  when  one  component  is  determined.  The  fourth  con¬ 
dition  is  provided  by  the  equilibrium  of  the  CO  water  shift 
reaction: 

CO  +  H20  -*  C02  +  H2  (18) 


Assuming  the  equilibrium  partial  pressure  should  be  P^o 

i  pe  pe 

H20’  rCCr  rCCh' 

-AGp\ 

RT  ) 


Pjjr  /jco>  /jco2-  we  have: 

Kp  =  P\\2  Pqq2  /  ^H,0  Px)  =  exP  ^ 


(19) 


where  A  Go  is  the  free  energy  change  of  the  shift  reaction 
when  all  reactants  and  products  are  at  P  =  1  atm.  After 
some  algebra,  one  finds: 

—b  +  sqrt  [b2  —  4(1  —  Kp)c] 


Pe  — 
rH2  - 


2(1  -  KP) 


(20) 


where 


b  =  Pco2  -  Pn2  +  Kp[PC0  +  2Ph2  +  Ph2o]  (21a) 

c  =  -[Ph2  +  Pu2o][Pco  +  Pu2]Kp  (21b) 


The  full  treatment  of  gas  flow  movements  requires  a  rig¬ 
orous  CFD  tool  and  is  computationally  demanding.  Though 
startup  and  transient  processes  as  well  as  variations  in  cer¬ 
tain  operating  parameters  may  have  a  sizeable  effect  on  the 
flow  profiles,  the  effect  on  the  overall  EC  performance  of 
the  cell  is  not  expected  to  be  of  the  same  order.  As  we  are 
primarily  concerned  with  the  thermal  and  EC  properties,  it 
is  not  necessary  to  solve  the  details  of  the  flow.  Instead,  it 
is  desirable  to  make  a  simplification  such  as  assuming  lam¬ 
inar  flow  to  reduce  the  computation  cost  and  allow  quick 
estimates  of  certain  flow  properties  that  are  influential  on 
the  thermal  and  EC  performance.  To  this  end,  the  most  rel¬ 
evant  data  are  the  species  concentration  profile  and  flow 
velocity. 

Given  the  current  density  profile  and  flow  rate,  the  species 
concentration  can  be  obtained  by  the  mass  balance  and  the 
boundary  conditions.  Assume  a  reaction  mixture  passes  at  a 
volume  rate  of  flow  u  (e.g.,  in  m3/s),  consider  an  element  of 
volume  dV  =  Adx  and  one  particular  species  component  k , 
which  enters  this  volume  element  at  a  concentration  Q,  and 
leaves  at  Ck  +  dC^  (in  mol/m3).  Here,  A  is  the  gas  channel 
cross-section  and  dx  is  the  distance  traveled  by  the  gas.  If 
there  is  no  longitudinal  mixing,  the  net  change  with  time  of 
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the  number  of  moles  of  k  within  dV,  dnk/dt ,  will  be  the  sum 
of  two  terms:  one  due  to  chemical  reaction  within  dV,  and 
the  other  equal  to  the  excess  of  species  k  entering  dV  over 
that  leaving.  Thus  we  have, 

d  nr 

-A  =  rkdV  -  udCk  (23) 

dt 


In  addition  to  the  concentration  and  pressure  changes  de¬ 
scribed  above,  the  pressure  drop  of  the  fuel  flow  through 
the  channel  due  to  viscosity  (assuming  laminar  flow)  is  es¬ 
timated  by: 


(R eDh) 


(30) 


The  chemical  reaction  rate  per  unit  volume  is  denoted  by  rk. 
For  H2  fuel,  we  may  write: 


rkdV  = 


Ids 

nF 


(24) 


where  a  is  the  fuel  channel  width  and  I  the  current  density. 
The  n  =  2  for  H2  and  n  =  4  for  O2  and  n  =  —2  for  H2O 
as  H2O  is  created  by  the  EC  reaction.  For  a  stable  species 
concentration,  dnk/dt  =  0,  we  have: 


la 


- dx  =  udCk 

nF 

(25) 

or 

fx  la 

Ck  —  Ck  0  /  dx 

Jo  unF 

(26) 

The  partial  pressure  can  be  calculated  as: 

Pk  =  CkRT=  Pk°  -  A  Pk 

(27) 

with 

fx  IaRT 

APk  =  /  v  dx 

Jo  unF 

(28) 

The  total  fuel  and  air  pressures  are  obtained  by  the  sum¬ 
mations  over  the  respective  species  partial  pressures.  For  H2 
fuel,  if  current  production  is  /,  the  rate  of  H2  mole  consump¬ 
tion  is  simply  dn\\2/dt  =  —I/(2F),  or 


dPH2  _  IRT 
dt  2  Fu 


(29a) 


where  Re  is  the  Dk -based  Reynolds  number,  Dk  the  hy¬ 
draulic  diameter,  /  the  length  of  the  flow  path,  and /depends 
on  the  shape  of  the  cross  section  of  the  channel,  e.g.,  /  = 
56.8  and  64  for  a  square  and  a  round  channel,  respectively 
[12].  Such  simplification  has  been  used  to  significantly  re¬ 
duce  the  computation  cost  [13]. 

2.5.  Heat  of  reaction 

If  one  is  interested  in  calculating  the  heat  of  chemical 
reaction  due  to  composition  change,  the  heat  of  reaction  is 
simply  the  difference  between  initial  system  heat  content 
(enthalpy)  and  the  final  system  heat  content.  If  one  is  con¬ 
cerned  with  EC  reaction,  then  the  heat  content  of  O2  needs 
to  be  included  in  the  calculation  of  heat  content  change. 
The  total  heat  production  is  the  heat  content  decrease  minus 
the  delivered  external  electrical  power,  Hi  —  Hf  —  V*pq, 
where  q  is  the  charge  involved  (product  of  current  and  time) 
and  Vop  the  fuel  cell  operating  voltage.  Hi  —  Hf  —  V*pq  is 
sufficient  for  a  2D  description,  which  is  one  option  of  the 
EC  module.  For  3D  models,  additional  details  are  the  fol¬ 
lowing:  Hi  —  Hf  —  q*  Vih  is  the  chemical  reaction  heat  and 
thermal  irreversibility  heat,  where  Vth  is  the  Nernst  poten¬ 
tial,  and  can  be  viewed  as  heat  released  in  the  fuel  channels 
and  anode.  q*(V th  —  VoP)  is  the  total  polarization  (Ohmic, 
transport,  activation)  heat  and  released  in  the  respective 
regions  in  accordance  with  polarization  loss.  In  other 
words,  the  overall  effective  “Ohmic”  resistance  of  the  cell 
is: 


For  H2  and  CO  fuel,  we  have: 


d  [Pr2  +  Pcq]  =  IRT 
dt  2  Fu 


(29b) 


The  shift  reaction  described  above  can  be  used  to  deter¬ 
mine  dP\\2/dt  and  dPco/dt.  For  fuel  with  CH4,  we  have: 

d[FH2  +  Pco  +  4FCh4]  _  IRT 

dt  ~  2  Fu  (  ° 

Together  with  the  equilibrium  composition  determination 
described  above,  the  respective  component  changes  can  be 
determined.  Once  the  species  concentrations  or  partial  pres¬ 
sures  are  known  at  the  gas  inlet,  the  above  equations  give 
the  variations  of  gas  compositions  and  partial  pressures  in 
the  gas  channels  when  the  temperature  and  current  density 
distribution  are  given. 

Notice  that  the  v  direction  is  the  flow  direction.  Other  than 
co-flow  configuration,  the  air  and  fuel  flows  have  different 
directions.  By  properly  accounting  for  the  flow  directions, 
all  co-,  counter-  and  cross-flow  designs  can  be  treated. 


Eoven  ~  V(i) 
R(i)  =  — — 

i 


=  Ri  + 


•  /?sinh 


-1 


In 


In 


^h2o*h2  / 


(31) 


The  total  heat  generated  is  i2R(i).  Heat  generated  by  each 
Ri  component  can  be  reasonably  assumed  to  be  released  to 
each  respective  component.  For  example,  heat  from  /o2-term 
is  released  in  cathode,  and  heat  generated  by  z'h2  and  z'h2o 
terms  are  released  in  the  anode.  Heat  from  b- term  is  assumed 
to  be  released  mainly  in  the  cathode  as  the  activation  polar¬ 
ization  is  dominated  by  the  cathode  process.  Therefore,  most 
activation  heat  should  be  released  in  cathode,  mainly  in  the 
reaction  zone  (about  10  p,m  thick)  at  the  cathode-electrolyte 
interface. 
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3.  The  EC  module  and  the  linkage  with  MARC® 

The  EC  module  for  the  fuel  cell  modeling  was  pro¬ 
grammed  in  a  stand-alone  code  using  FORTRAN  77  and 
can  be  executed  if  temperature  field  and  the  inlet  flow  in¬ 
formation  are  available.  The  module  is  arranged  as  a  single 
user-defined  function  (UDF)  for  interaction  with  external 
software.  The  main  program  structure  can  be  seen  as  the 
following: 


The  UDF  also  computes  the  current  density  for  a  specified 
voltage.  The  current  density  is  obtained  by  bisection  search 
method.  The  initial  current  density  range  is  between  zero 
and  the  smaller  of  cathode  and  anode  limiting  current  (less 
a  small  quantity,  10-6,  to  avoid  the  logarithmic  singularity). 
The  implementation  of  several  parameters  and  two  forms 
of  I-V  curves  allows  for  the  code  to  be  calibrated  easily  to 
available  data  taken  from  SOFC  testing.  The  user  can  also 
create  other  desired  I-V  formula. 


The  EC  code  functions  in  several  modes.  These  include: 

•  For  a  given  operating  voltage  of  a  cell  plate,  the  model 
computes  the  current  density  distribution  according  to  the 
given  uniform  cell  voltage. 

•  For  a  given  total  current,  the  model  finds  the  cell  oper¬ 
ating  voltage  so  that  the  integration  of  the  local  current 
density  with  the  operating  voltage  yields  the  required  total 
current. 

•  For  a  stack  operation  and  given  total  voltage,  the  model 
determines  each  individual  cell  voltage  so  that  they  sum  to 
the  desired  total  voltage  and  the  total  current  is  the  same 
for  each  individual  cell. 

3.1.  The  I-V  sub-model 

One  function  of  the  UDF  is  to  calculate  current  density, 
cell  voltage,  and  heat  production  in  an  SOFC  stack,  tak¬ 
ing  as  inputs  local  values  of  the  gas  partial  pressures  and 
temperatures.  The  details  of  I-V  relation  and  parameters 
therein  have  been  described  above.  For  recap,  the  cell  po¬ 
tential,  V(7),  for  a  given  current  density,  given  also  the  local 
values  of  the  temperatures  and  gas  partial  pressures,  is  de¬ 
termined  by  the  chosen  I-V  relation.  It  has  the  following 
characteristics: 

(a)  Temperature  and  geometry  dependence  of  Ohmic  elec¬ 
trical  resistance  is  considered  explicitly.  One  is  free  to 
adjust  the  thickness  of  each  individual  cell  component. 

(b)  Temperature  dependences  of  non-Ohmic  overpotential 
terms  by  interpolations,  aided  by  a  few  experimental 
data. 

(c)  Partial  pressure  and  geometry  dependences  of  the 
non-Ohmic  terms  are  included  according  to  microscopic 
theory. 


3.2.  Grid 

To  adequately  describe  the  local  distribution  of  current 
density,  heat  generation,  and  species  concentrations,  a  3D 
grid  is  generated  for  a  fuel  cell(s)  in  a  stack: 

(a)  In  the  plan  of  the  cell  plate,  a  rectangular  mesh  is  used 
for  the  EC  active  area.  The  grid  points  (v,  y),  across  the 
PEN  structure,  z,  are  associated  pair-wise.  Tests  show 
one  grid  point  for  each  (1  cm  x  1  cm)  area  gives  adequate 
EC  results,  while  one  point  for  each  (0.5  cm  x  0.5  cm) 
area  provides  practically  convergent  results. 

(b)  For  a  set  (x,  y)  in  a  cell,  current  is  constant  over  the 
z-grid.  The  I-V  relation  is  used  to  compute  cell  voltage 
if  total  current  density  is  given  or  to  find  the  current 
density  consistent  with  a  given  cell  voltage. 

(c)  All  cells  in  the  stack  have  the  same  grid,  but  correspond¬ 
ing  grid  points  in  different  cells  can  have  different  cur¬ 
rent  densities. 

(d)  Distributed  heat  productions  are  determined  over  the 
grids  for  the  various  terms  in  the  I-V  relation  according 
to  their  respective  physical  origins. 

(e)  As  fuel  cells  are  much  thinner  than  the  cell  plate  dimen¬ 
sion,  very  few  z  points  are  required  to  obtain  satisfying 
results.  In  fact,  one  can  also  choose  to  use  a  2D  grid, 
i.e.,  one  z  point  for  each  PEN  structure. 

3.3.  Cell  voltage  for  a  required  total  current 

When  the  total  current  through  a  cell  has  been  specified, 
the  potential  across  the  cell  is  determined  according  to  the 
equipotential  condition  and  resulting  current  density  distri¬ 
bution  integrated  to  the  specified  total  current.  The  search  is 
done  by  bisection  with  an  initial  range  of  between  zero  and 
the  open  circuit  voltage. 
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(a)  The  search  for  current  for  a  voltage  in  the  I-V  submodel 
described  above  is  used  as  a  subset.  The  search  is  done 
for  every  (v,  y)  node. 

(b)  Current  at  each  node  is  integrated  and  compared  to  the 
specified  total  current.  The  search  continues  until  con¬ 
vergence  is  achieved. 

(c)  The  search  for  each  (v,  y)  node  is  independent  of  the 
others,  and  the  total  cost  is  linear  to  the  number  of  (v,  y) 
nodes.  Cost  of  integration  (step  b)  above  is  also  linear 
to  the  number  of  (v,  y)  nodes.  Therefore,  the  total  cost 
is  linear  to  the  number  of  (x,  y )  nodes. 

3.4.  Self-consistence  for  the  current  generation  and 
concentration  distribution 

The  existing  I-V  subroutines  require  temperature  profile 
and  gas  pressure  distributions  in  the  channels  as  input  to  find 
the  desired  I  or  V.  In  practice,  only  gas  partial  pressures  in 
the  gas  inlets  are  known  working  parameters.  The  variation 
of  partial  pressures  inside  the  channels  is  determined  by  the 
current  density  distribution,  as  shown  in  the  subsection  of 
the  flow  model  (Fig.  1).  Therefore,  there  is  a  coupled  pro¬ 
cess  between  the  current  distribution  and  gas  pressure  distri¬ 
butions.  A  self-consistent  determination  of  the  gas  pressures 
is  required.  This  is  done  in  the  code  by: 

(a)  An  initial  guess  of  the  current  density  distribution  is 
generated  by  the  overall  operating  conditions  (e.g.,  inlet 
fuel  pressures,  working  temperatures,  cell  voltage,  or 
total  current). 

(b)  The  current  density  distribution  is  used  to  find  the  gas 
pressure  changes  in  the  channels,  as  shown  in  the  flow 
model. 

(c)  The  resulting  gas  pressures  are  used  to  determine  the 
new  current  density. 


(d)  The  new  current  density  distribution  is  mixed  with  the 
old  current  density  to  get  the  revised  guess  current  den¬ 
sity  distribution. 

(e)  Compare  the  total  current,  if  operating  voltage  is  set, 
or  the  operating  voltage,  if  the  total  current  is  speci¬ 
fied,  and  iterating  until  the  self-consistent  convergence 
is  achieved. 

3.5.  Stack  operation 

(a)  When  the  total  current  through  a  stack  of  cells  is  spec¬ 
ified,  individual  cell  potential  is  determined  separately 
and  their  sum  gives  the  total  output  voltage.  Each  in¬ 
dividual  cell  is  independent  of  the  others.  The  CPU  re¬ 
quirement  is  linear  to  the  number  of  cells  in  the  stack. 

(b)  If  the  total  voltage  is  to  be  specified  for  a  stack,  all 
the  cell  voltages  are  adjusted  such  that  their  sum  is  the 
specified  value  while  the  total  current  is  the  same  for 
each  cell.  As  total  current  should  be  invariant  for  all  the 
cells,  it  is  used  as  the  searching  variable  (i.e.,  case  (a)  is 
used  as  a  subset).  The  cost  is  again  linear  to  the  number 
of  cells  in  the  stack.  Typically,  case  (b)  requires  10-20 
times  more  CPU  time  than  case  (a). 

The  above  is  done  self-consistently  and  the  resulting  cur¬ 
rent  distribution,  heat  productions  and  cell  voltages  are  all 
returned  to  the  user. 

3.6.  EC  module  summary 

We  have  described  the  theory  and  implementation  details 
of  an  EC  model  capable  of  computing  the  current  density 
distribution,  chemical  and  electrical  heat  generation,  and 
fuel  cell  efficiency.  The  EC  model  can  address  the  depen¬ 
dence  of  fuel  cell  performance  on  the  operating  tempera- 


Fig.  1.  MARC®  model  showing  different  materials  through  the  thickness. 
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ture,  fuel  composition,  gas  pressure,  flow  rate,  and  geome¬ 
try/configuration.  The  module  can  be  linked  easily  to  other 
simulation  software  to  improve  modeling  quality,  e.g.,  pro¬ 
viding  the  FEA  code  with  EC  heat  flux  to  facilitate  ther¬ 
mal  stress  analysis,  or  assisting  CFD  with  heat  term  and  gas 
composition  change. 

3.7.  Linking  the  EC  module  to  MARC® 

The  basic  procedure  for  the  FEA-EC  integration  was 
to  develop  intermediary  subroutines  that  will  pass  the  cell 
temperature  distribution  to  the  EC  subroutine  and  return 
the  heat  generation  to  MARC®.  This  is  accomplished  with 
a  FLUX  subroutine  and  it  is  supported  by  the  use  of  state 
variables.  State  variables  are  used  to  track  the  composition 
of  the  gases;  H2,  H2O,  CO,  CO2  and  CH4  in  the  fuel,  and 
N2  and  O2  in  the  air.  Additional  state  variables  track  the 
heat  generation  and  current  density,  which  can  be  plotted 
subsequently  through  the  plotting  subroutine.  The  interface 
can  be  viewed  simply  as 

MARC  source  subroutine  (flux(f,ts,n,time)) 


that  the  heat  conduction  is  equivalent  to  the  actual  convective 
heat  transfer.  The  anode,  electrolyte,  and  cathode  are  each 
modeled  with  a  single  element.  Elements  representing  the 
interconnect  are  located  at  the  top  and  bottom  of  the  stack. 
Some  of  the  material  properties  and  cell  component  dimen¬ 
sions  are  shown  in  Table  1.  Other  parameters  include  the 
film  coefficient:  hair  =  102.3  J/(sK);  hH2  =  1530J/(sK); 
corresponding  to  an  equivalent  conductivity  of  film  layer  = 
h (film  coefficient)  x  £(film  layer  thickness).  Heat  transfer 
mechanisms  considered  are  convection  and  conduction.  The 
model  uses  a  10  x  10  mesh  of  elements  to  represent  the 
in-plane  dimensions  of  the  cell,  which  is  an  11  cm  x  1 1  cm 
square  cell.  For  a  multi-cell  stack,  one  can  build  a  similar 
procedure  file  by  repeating  the  corresponding  layers. 

To  illustrate  the  modeling  capability  of  the  MARC®-EC 
package,  Fig.  2  shows  the  steady-state  temperature  profiles 
for  three  flow  designs  (co-,  counter-  and  cross-flow)  with 
H2  +  CO  fuel.  Comparing  the  results  with  two  different  fuel 
flow  rate  (Vfuei  =  0.055  and  0.0275 1/s,  respectively),  one 


Input  temperature  profile  and  geometry  info  4  freturn  heat  flux  and  other  state  variables 

EC  module  interface  subroutine 


The  EC  subroutine  is  self-consistent  and  solves  the  en¬ 
tire  mesh  in  one  pass.  The  FLUX  subroutine  is  called  at 
every  integration  point  for  every  element  in  the  PEN.  This 
discrepancy  was  accommodated  by  executing  the  EC  sub¬ 
routine  from  FLUX  only  once  at  the  first  integration  point 
of  the  first  element  and  placing  the  results  in  an  array.  Sub¬ 
sequent  calls  to  FLUX  then  simply  read  the  appropriate 
value  stored  in  the  array. 

To  use  the  MARC®-EC  package,  one  is  required  to  pro¬ 
vide  the  necessary  information  such  as  the  SOFC  geometry, 
material  properties,  boundary  conditions,  operation  param¬ 
eters  (e.g.,  fuel-air  flow,  working  temperature,  etc.).  Such 
information  can  in  principle  be  given  through  the  MARC® 
graphical  interface.  It  can  also  be  specified  in  a  procedure 
file.  These  procedures  are  common  in  using  MARC®,  and 
readers  are  referred  to  the  MARC®  documentation  for 
details. 


4.  Sample  results 

As  an  example,  a  simple  one-cell  model  is  shown  in 
Fig.  1 .  The  model  uses  1 1  elements  through  the  thickness  of 
the  cell  with  the  following  1 1  layers:  two  interconnect  layers 
(SS430),  two  fuel  film  layers,  two  air  film  layers,  one  fuel 
layer,  one  air  layer,  one  cathode  layer,  one  electrolyte  layer 
and  one  anode  layer.  The  fuel  and  air  are  represented  by 
three  elements  each.  Thinner  elements,  designated  as  film  el¬ 
ements,  are  located  between  the  fluid  and  the  solid  elements. 
The  thermal  conductivity  of  these  elements  is  adjusted  so 


can  see  the  significant  impact  of  flow  rate  on  the  temperature 
distribution  (and,  therefore,  on  the  thermal  stress). 

The  modeling  capability  of  the  MARC®-EC  tool  is  not 
limited  to  the  steady-state  analysis.  This  modeling  tool  can 
also  be  used  for  the  evaluation  of  the  transition  from  startup 
to  steady-state  conditions.  As  an  example,  Fig.  3  shows  some 
snapshots  of  a  fuel  cell  from  startup  to  steady- state  transi¬ 
tion.  The  fuel  cell  was  initially  at  T  =  0  °C.  The  fuel  and  air 
flow  were  at  T  =  700  °C  and  with  flow  rate  of,  respectively, 
0.0825  and  0.33 1/s.  The  fuel  cell  was  set  at  a  working  volt¬ 
age  of  V  =  0.7  V  when  it  is  warm  enough  to  have  current 
production. 

Due  to  the  nature  of  the  MARC®  FEA  method,  the  numer¬ 
ical  results  are  very  stable.  The  entire  comparable  transient 
results  are  almost  identical  for  time  steps  ranging  from  0.2  to 
12  s — a  significant  advantage  over  a  typical  CFD  algorithm. 
The  startup  modeling  has  demonstrated  that  startup  will  not 
occur  in  an  isothermal  manner.  Steady- state  operating  con¬ 
ditions  will  result  in  a  different  temperature  distribution  due 
to  the  heat  transfer  and  fuel  consumption  along  the  length 
of  the  cell.  The  startup  procedure  and  steady-state  condi¬ 
tions  will  each  need  be  controlled  to  maintain  an  adequate 
safety  factor  below  the  material  failure  points.  In  the  same 
manner,  however,  starting  the  EC  reaction  after  the  cell  is 
heated  also  needs  to  be  carefully  controlled.  This  tool  can 
be  used  to  identify  and  prevent  conditions  where  localized 
thermal  runaway  may  occur.  Similarly,  any  disturbance  to 
steady- state  operating  conditions  can  be  evaluated. 

As  another  example,  Fig.  4  shows  the  cell  voltage  vari¬ 
ation  for  a  30-cell  stack  (Vtot  =  21V)  at  the  steady  state. 
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Fig.  2.  Temperature  profiles  of  cross-,  co-  and  counter-flow  designs.  Upper  row:  Vfuei  =  0.055 1/s;  lower  row:  Vfuei  =  0.0275 1/s 
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Fig.  3.  Snapshot  of  the  transition  from  start-up  to  steady-state  conditions.  Shown  are  for  start-up  time,  t :  (a)  12s,  (b)  60s,  (c)  300s,  (d)  600s,  (e)  900s 
(f)  1200  s. 
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Table  1 


Some  material  and  geometry  parameters  used  in  the  simulation 


Thickness  (m) 

Heat  conductivity  (J/mKs) 

Specific  heat  capacity  (J/kgK) 

Mass  density  (kg/m3) 

Interconnect  layer 

5e  —  4 

25.5 

800 

7700 

Air  film  layer 

le  —  4 

102.3e  -  4 

100 

0.3236 

Air  layer 

1.5e  —  3 

0.0692 

1156 

0.3236 

Cathode 

5e  —  5 

2.5 

520 

6350 

Electrolyte 

le  —  5 

2 

460 

6010 

Anode 

6e  —  4 

22.5 

466 

4200 

Fuel  film  layer 

le  —  4 

1530e  -  4 

100 

0.0224 

Fuel  layer 

5e  —  4 

0.450 

15156 

0.0224 

cell  voltages  for  a  30  cell  stack  (Vtot=21V) 


cell  number 

Fig.  4.  The  cell  voltage  variation  of  a  30-cell  stack  (Vtot  =  21  V)  at  the  steady  state. 


As  can  be  seen,  only  the  outmost  top  and  bottom  2-3  cells 
deviate  substantially  from  the  average  cell  voltage.  Similar 
results  are  found  for  8-  and  15-cell  stacks. 

Numerical  characteristics  of  the  MARC®-EC  tool  can  be 
seen  in  Tables  2  and  3  (performed  on  one  SGI  RISC  14000 
processor,  400  MHz).  The  CPU  time  requirement  increases 
roughly  linearly  with  the  number  of  cells  in  the  stacks.  Typ¬ 
ically  20-40  increments  are  needed  for  a  steady-state  run. 
Therefore,  one  can  get  the  steady-state  results  for  a  30-cell 
stack  in  a  few  minutes  if  the  total  current  is  the  desired  cri¬ 
teria.  Even  if  the  total  stack  voltage  is  the  control  criteria, 
the  results  can  be  obtained  in  a  few  hours  of  CPU  time  in  a 

Table  2 


CPU  time  when  specifying  total  current  or  fuel  utilization 


No.  of  cells 

Time/inc  (s) 

2 

0.45 

3 

0.68 

4 

0.91 

5 

1.15 

8 

1.83 

15 

3.48 

30 

6.93 

Table  3 


CPU  time  when  specifying  total  voltage 


No.  of  cells 

Time/inc  (s) 

2 

14.0 

3 

26.8 

4 

36.3 

5 

51.8 

8 

53.0 

15 

149 

30 

326 

single  process.  Consequently,  one  can  examine  various  stack 
designs  readily  with  this  tool. 

5.  Summary 

We  have  described  the  technical  details  of  a  multi-physics 
computational  tool  for  SOFC  modeling,  the  MARC®-EC 
package.  The  description  focuses  on  the  development  of 
an  electrochemistry  module  to  supplement  the  capability 
of  MARC®  to  model  SOFCs.  Numerical  results  are  given 
to  demonstrate  the  modeling  capability  and  numerical 
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efficiency  of  the  tool.  The  tool  can  be  used  to  perform 
parametric  studies  that  are  important  to  system  designs. 
The  MARC®-EC  package  offers  significant  advantages  in 
terms  of  numerical  stability  and  computational  efficiency 
over  that  of  a  CFD  code  such  as  STAR-CD  or  Fluent. 
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